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Tracking  performance  depends  upon  the  quality  of  the  measurement  data. 

In  the  Kalman-Bucy  filter  and  other  trackers,  this  dependence  is  well- 
understood  in  terms  of  the  measurement  noise  covariance  matrix,  which 
specifies  the  uncertainty  in  the  values  of  the  measurement  inputs.  When 
the  origin  of  the  measurements  is  also  uncertain,  one  has  the  widely- 
studied  problem  of  data  association  (or  data  correlation) ,  and  tracking 
performance  depends  critically  on  additional  parameters,  primarily  the  — 
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ABSTRACT  (continued) 


probabilities  of  detection  and  false  alarm.  In  this  paper  we  derive 
a  modified  Riccati  equation  that  quantifies  (approximately)  the 
dependence  of  the  state  error  covariance  on  these  parameters.  We 
also  show  how  to  use  an  ROC  curve  in  conjunction  with  the  above 
relationship  to  determine  an  "optimal"  detection  threshold  in  the  signal 
processing  system  that  provides  measurements  to  the  traclcer. 
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1 .  INTRODUCTION 


Garden-variety  tracking  problems  involve  processing 
measurements  (e.g.,  range  and  azimuth  observed  by  a  sensor)  from 
a  target  of  interest  and  producing,  at  each  time  step,  an 
estimate  of  the  target^’s  current  position  and  velocity  vectors. 
Uncertainties  in  the  target  motion  and  in  the  measured  values, 
usually  characterized  as  random  noise,  lead  to  corresponding 
uncertainties  in  the  target  state. 

A  common  and  versatile  approach  to  such  problems  involves 
assuming  that  the  state  dynamics  and  the  measurements  are  both 
corrupted  by  additive,  white,  possibly  Gaussian  noise;  the 
solution  is  then  the  celebrated  Kalman-Bucy  filter 

[1,  2,  3,  4,  5],  which  is  the  conditional  mean  state  estimator, 
best  linear  estimator,  maximum  a  posteriori  estimator,  maximum 
likelihood  estimator,  or  least-squares  estimator,^  depending  upon 
oner's  point  of  view.  The  parameters  that  determine  tracking 
performance  in  such  a  filter  are  the  system  matrices  in  the 
equations  describing  target  state  dynamics  and  measurements, 
which  will  be  considered  fixed  for  the  purposes  of  this 
discussion,  and  the  covariance  matrices  of  the  process  and 
measurement  noise,  which  specify  the  uncertainties  in  target 
motion  and  measured  values,  respectively. 

In  many  tracking  problems,  particularly  those  arising  in 
surveillance,  there  is  additional  uncertainty  regarding  the 
origin  of  the  received  data,  which  may  (or  may  not)  include 


^In  the  least-squares  case,  the  assumptions  about  noise  are 
replaced  by  assumptions  about  error  weightings. 


measurements  from  the  target (s)  of  interest,  interfering  targets, 
or  random  clutter  (false  alarms) .  This  leads  to  the  problem  of 
data  association  or  data  correlation,  which  has  been  attacked  on 
a  number  of  fronts  [6,  7,  8,  9,  10,  11,  12,  13,  14]  and  surveyed 
in  [15,  16,  17].  In  this  situation,  tracking  performance  depends 
not  only  upon  the  noise  covariances,  but  upon  the  amount  of 
uncertainty  in  measurement  origin.  In  some  of  the  approaches 
cited  above  [6,  7,  8,  9,  10],  this  dependence  is  explicit  and  is 
characterized  in  terms  of  the  detection  probability  and  clutter 
density  (which  is  proportional  to  false  alarm  probability) . 

In  typical  applications,  measurement  data  are  provided  to  a 
tracker  by  upstream  signal  processing  and  detection  algorithms, 
as  indicated  in  figure  1.  The  process  noise  covariances  are 
normally  selected  on  the  basis  of  experience  and  intuition  (i.e., 
they  are  guessed) .  The  measurement  noise  covariances  are  either 
provided  by  the  signal  processing  algorithm,  as  shown  in  the 
figure,  or  they  are  selected  in  the  same  manner  as  the  process 
noise.  In  any  case,  the  true  noise  levels  are  usually  fixed  by 
target  dynamics  and  sensor  configuration  and  cannot  be  adjusted 
on  line. 

Detection  and  false  alarm  probabilities,  on  the  other  hand, 
are  highly  interdependent  and  adjustable  via  a  de tection 
threshold;  raising  the  threshold  lowers  both  probabilities,  and 
vice-versa.  This  relationship,  which  is  also  parameterized  by 
signal-to-noise  ratio  (SNR) ,  is  usually  characterized  by  means  of 
a  set  of  receiver  operating  characteristic  (ROC)  curves,  as 
discussed  below  in  Section  4.  The  threshold  is  typically  set  by 
choosing  a  design  point  on  the  most  applicable  ROC  curve,  based 
on  the  perceived  tradeoffs  between  false  alarms  and  missed 
detections.  However,  to  the  best  of  our  knowledge,  these 
tradeoffs  have  never  included  any  systematic  or  quantitative 
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FIG.  1.  TRACKING  SYSTEM  BLOCK  DIAGRAM 


consideration  of  the  effects  downstream  on  data  association  and 
tracking  performance. 

In  this  paper  we  shall  derive  such  a  quantitative 
relationship,  in  which  the  dependence  of  a  tracker^'s  error 
covariance  upon  detection  and  false  alarm  probability  is 
explicitly  (but  approximately)  characterized  by  a  scalar 
parameter  in  the  covariance  equation  (modified  Riccati  equation) . 
This  result  is  important  for  the  following  reasons: 

1.  Contour  plots  of  the  scalar  parameter  as  a  function  of 
detection  probability  and  false  alarm  probability  form 
a  set  of  "tracker  operating  characteristic"  (TOC) 
curves which  can  be  superimposed  on  ROC  curves  for  the 
detector  or  receiver  of  interest  in  order  to  determine 
graphically  the  operating  point  that  "optimizes" 
tracker  performance. 

2.  The  stability  of  the  tracking  process  depends 
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critically  on  the  detection  and  false  alarm 
probabilities;  whether  the  modified  Riccati  equation 
remains  stable  for  all  values  of  the  scalar  parameter 
remains  an  open  question. 

3.  Allocation  of  tracking  resources  (both  computation  and 
communication)  requires  prediction  of  future  state 
error  covariances  under  various  resource 
configurations,  i.e.,  as  a  function  of  detection  and 
false  alarm  probability  and  of  process  and  measurement 
noise  covariance. 

4.  The  same  derivations  provide  a  solution  to  the  related 
problem  of  determining  the  statistical  properties  of 
the  modified  likelihood  function  [18] ,  used  for 
decision  making  (e.g.  maneuver  detection)  when 
measurement  origins  are  uncertain. 

In  Section  2  the  problem  of  relating  tracking  performance  to 
detection  and  false  alarm  probabilities  is  formulated  in  the 
context  of  probabilistic  data  association.  The  key  results  are 
derived  in  the  next  section,  followed  by  analyses  and  ROC  curves 
of  simple  receivers  in  Section  4.  These  results  are  combined  in 
Section  5  and  used  to  select  optimal  operating  points  on  the  ROC 
curves.  Application  to  the  modified  likelihood  function  is 
discussed  in  the  next  section,  followed  by  concluding  remarks  in 
Section  7. 
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2 .  PROBLEM  FORMULATION 


Consider  a  dynamic  system  (target  model)  of  the  usual  form, 
ik+i  =  £Xk  +  (1) 

Zk  =  Mk  +  ^k  .  (2) 

where  x  is  the  state  vector,  ^  is  the  measurement  vector,  w  and  v 
are  zero-mean,  mutually  independent,  white,  gaussian  noise 
vectors  with  covarian.f"  matrices  Q  and  R,  respectively,  and  k  is 
a  discrete  time  index.  The  matrices  F,  G,  H,  Q,  and  R  are 
assumed  known  and  their  dependence  on  t  is  suppressed  here  for 
notational  convenience.  The  initial  state  is  assumed  gaussian 
with  mean  ^o|0  covariance  Po|0*  ^  typical  state  vector  would 

include  position  and  velocity  variables,  as  well  as  other 
information  that  relates  to  the  specific  type  of  platform  being 
tracked,  and  typical  dynamics  would  assume  straight-line  or 
great-circle  target  motion  with  disturbances  from  the  process 
noise  w. 

The  tracker's  estimate  of  the  target  state  Xy.  at  time  k, 
given  data  up  to  time  i,  is  denoted  ^kji*  error  in  this 

estimate  is  =  2ik”^|i'  with  error  covariance  matrix  Pk  |  i  = 

k|i^'  where  E  denotes  expectation.  The  discrete-time 
Kalman-Bucy  filter  [2,  3,  4,  5]  propagates  these  in  two  stages. 
The  prediction  stage  accounts  for  time  evolution, 

ikjk-l  “  l!ik-llk-l  ^2) 

Zklk-l  “  ££k-l|k-ir  +  GQG'  (4) 
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starting  from  the  initial  conditions  1 0  Zo|0‘  update 

stage  compares  the  incoming  measurement  with  the  predicted 
measurement  2j^  *  1 1-1  form  the  innovation  vector 

2k  -  Zk  “  £k|k-l  (5) 

whose  covariance  matrix  is 

Sk  =  E{2k2k}  =  HPklk-ir  +  R  (6) 

The  state  and  covariance  are  then  updated  via 

^Iclk  *  ik|k-l  Skik 

£k|k  *  <I-Wka)£k|k-l(i-«k2)'  +  «k2«k 
=  £k|k-l  -  WkSk!!k 

where  Wj^  =  2k|k-  is  the  gain  matrix.  The  resulting  state 

estimate,  under  the  above  assumptions,  is  the  conditional  mean 
— k  I  k  “  e{x|^|Y*^}  where  denotes  all  data  vectors  for  i  £  k. 

In  order  to  avoid  clouding  the  discussion  to  follow,  this 
brief  summary  ignores  a  number  of  complications  that  arise  in 
practice.  If  the  system  is  nonlinear,  for  example,  then  it  can 
usually  be  linearized  and  the  same  basic  equations  can  be  applied 
to  deviations  from  the  nominal  trajectory  [3,  5].  If  the  target 
occasionally  deviates  from  the  assumed  motion  model,  e.g.,  by 
maneuvering,  then  some  decision-making  or  other  machinery  must  be 
provided  to  deal  with  these  instances. 

In  multi-sensor  problems,  the  size  and  composition  of  the 
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measurement  vector  often  varies  from  one  time  to  the  next;  in 
other  words,  is  composed  of  independent  subvectors  from 
various  sensors,  any  subset  of  which  may  be  present  at  a  given 
time.  Moreover,  in  the  problem  of  interest  here,  each  sensor 
supplies  not  one  but  several  subvectors  that  must  be  associated 
with  targets.  We  will  avoid  the  resulting  notational  morass  by 
restricting  equations  (5) -(8)  to  apply  to  a  measurement  subvector 
Xjt  from  a  single  sensor.  In  addition,  we  will  suppress  the  time 
subscript  k  from  all  variables  except  P  and  Y,  unless  it  is 
required  for  clarity.  Without  any  loss  of  generality,  the  data 
association  problem  may  now  be  formulated  as  follows. 

At  each  time  step,  the  sensor  provides  a  set  of  candidate 
measurements  to  be  associated  with  targets  (or  rejected) .  In 
most  approaches,  this  is  done  by  forming  a  "validation  gate" 
around  the  predicted  measurement  from  each  target  and  selecting 
those  detections  that  lie  within  the  gate.  There  are  many 
different  approaches  to  establishing  a  correspondence  between 
candidate  measurements  and  targets;  some  of  these  were  cited 
above  in  Section  1. 

In  this  paper  we  shall  focus  on  the  probabilistic  data 
association  (PDA)  method  [6,  7,  15],  although  the  results  are 
relevant  to  other  methods  [9,  10]  in  which  similar  machinery  is 
used.  The  m  candidate  measurements  are  denoted  ^ j ,  j=l,...m,  and 
their  corresponding  innovations  are 

2j  “  Zj  -  If  j=l,...m  (9) 

The  term  measurement  will  be  used  interchangeably  for  and  2j , 
since  they  contain  equivalent  information  [5] . 

Considering  a  single  target  independently  of  any  others,  Xj 
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denotes  the  event  that  the  j-th  measurement  belongs  to  that 
target  and  Xg  the  event  that  none  of  the  measurements  belongs  to 
it  (no  detection) .  The  PDA  approach  builds  upon  the  assumptions 

that  the  estimation  errors  x  and  £  have  gaussian  densities  at 
each  time  step  (this  is  approximate,  since  there  is  an 
exponentially  growing  tree  of  possible  measurement  sequence 
hypotheses  and  the  true  densities  are  weighted  sums  of 
gaussians) .  It  is  also  assumed  that  the  correct  measurement  is 
detected  with  probability  (independently  at  each  time)  and 

that  all  other  measurements  are  Poisson-distr ibuted^  with 
parameter  CV,  where  C  is  the  expected  number  of  false 
measurements  per  unit  volume  and  V  is  the  volume  of  the 
validation  gate. 

The  gate  is  normally  a  g-sigma  ellipsoid  {£  ;  )  and 

Pq  is  the  probability  that  the  correct  measurement,  if  detected, 
lies  within  the  gate.^  The  gate  volume  is  thus  V  « 
where  M  is  the  dimension  of  £  and  Cj^=TrM/2/p (jj/2+1)  j-he  volume 

of  the  M-dimensional  unit  sphere  (Cj^=2,  C2*tt,  C3=4’nr/3,  etc.). 

The  conditional  mean  estimate  ^  is  obtained  from  (7)  by 
using  the  combined  (weighted)  innovation 


(10) 


^Equivalently,  the  number  n  of  false  measurements  has 
probability  mass  function  p(n)  =  e"^^(CV)”/nl  and  the  location  of 
each  false  measurement  is  uniformly  distributed  in  the  gate. 

%his  is  just  the  gaussian  probability  mass  in  the  gate,  which 
is  pfien  assumed  to  be  unity  in  practice,  since  Pg>.99  whenever 
q>Vl^'^+2. 
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where  8j  =  PiXjlY*'},  is  the  posterior  probability 

that  the  j-th  measurement  (or  no  measurement,  for  j*0)  is  the 
correct  one.  These  probabilities  are  given  by  the  following 
expressions  (see  Appendix  A) : 


exp(-£jS~^£j/2) 
b  +  ^  exp(-2^s'^£i/2) 


b 


b  + 


exp(-£^S"^£i/2) 


3“1, • • .m 


(11) 


(12) 


where 

b  S  (2Tr)M/2c|s|  V2(i_p^p^)/p^ 

=  (2Tr)M/2(CV/c„gM)  (1-Pj3Pq)/Pj3  (13) 


The  covariance  update  equation  (8)  becomes 

ik|k  =  Pk|k-1  -  (i-8o)wsw'  +  Zk  (14) 

where 

Zk  “  Wt2^Sj£j£j-££nr  (15) 

The  data-dependent  factors  6q  and  P|^  transform  the  original 
deterministic  Riccati  equation  into  a  stochastic  one. 
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The  principal  purpose  of  this  paper  is  to  explore  the 
dependence  of  tracking  performance,  particularly  the  behavior  of 
on  detection  probability  Pjj  and  clutter  density  C.  This  is 
accomplished  in  the  next  section  by  means  of  a  deterministic 
approximation  to  the  stochastic  nidified  Riccati  equation  (14) . 

Another  major  problem  in  data  association  and  tracking  is 
the  testing  of  hypotheses  for  maneuver  detection,  track 
initiation,  signature  formation,  target  classification,  and  other 
decision-making  purposes.  The  uncertainty  in  measurement  origin 
leads,  in  the  PDA  and  other  approaches,  to  a  modified  likelihood 
function  [18]  involving  the  combined  innovation  2.  ^ 

major  drawback  of  this  approach  was  the  difficulty  of  computing 
the  covariance  matrix  of  £»  this  can  now  be  done  using  an 

intermediate  result  to  be  derived  in  the  next  section.  Its 
application  to  the  modified  likelihood  function  will  be  discussed 
in  Section  6. 

Finally,  note  that  multiple  targets  can  be  handled 
simultaneously  via  the  joint  probabilistic  data  association 
(JPDA)  approach  (8] ,  in  which  the  posterior  probabilities 
(11) -(12)  are  computed  jointly  across  a  set  of  potentially 
interfering  targets.  Although  this  is  a  very  important 

extension,  it  greatly  complicates  the  derivations  and  will  not  be 
included  here. 
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3.  APPROXIMATE  COVARIANCE  EQUATION 

Since  any  measure  o£  tracking  performance  must  depend 
heavily  (or  perhaps  exclusively)  on  the  error  covariance  matrix 
P|^|l^,  we  shall  attempt  to  characterize  its  behavior  in  the 
presence  of  uncertainties  in  measurement  origin.  is  a 

random  (data-dependent)  matrix  governed  by  the  nonlinear, 
stochastic  difference  equations  (4)  and  (14)  ,  and  hence  its 
behavior  can  only  be  determined  in  a  statistical  sense. 
Moreover,  even  propagation  of  its  first  and  second  moments 
appears  to  be  intractable  except  via  extensive  numerical 
operations . 

Consequently,  we  shall  consider  an  approximation  to  (14)  in 
which  the  random  matrix  ^  defined  in  (15)  and  the  probability  6q 
given  by  (12)  are  replaced  by  their  (prior)  expected  values 

Zk  »  eIPj^Iy''”^}  (16) 

Bo  -  e{boIy''"^}  -  e{6o}  •  1  -  PdPg  (1-7) 

where  the  final  expression  is  a  consequence  of  E [p{a| B} ] «p{a} . 

These  substitutions  make  (4)  and  (14)  into  a  set  of 
deterministic  equations  that  can  be  iterated  forward  in  time. 
Because  (14)  is  nonlinear  in  P|5|)5_i,  this  does  not  yield  e{P]^||^}} 
nevertheless,  it  will  give  approximate  values  of  future  state 
error  covariances  in  the  presence  of  uncertain  detections  and 
false  alarms  as  a  function  of  the  environmental  parameters  Pq  and 
C,  and  of  the  noise  covariances  R  and 

Expansion  of  (16)  yields 
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Ik  -  -  ElElP^lm^Y*^-^] 

00 

m»0 

where  p{iii|Y**”^}  »  P{«}  is  given  by  (37)  in  Appendix  A.  Using 
(15)  and  (10) ,  the  inside  expectation  becomes 


E(P|jlm,Y'^-’^] 


W(U]^(m)-U2(m)  IW',  m»l,2, 


0, 


m«0 


where 


Ui(m)  a  E(^  Siiiiilm^Y**"!) 

I 

U2(m)  ^  ®i2i^  ®j2jl>nrY>'-^] 


(19) 


(20) 


(21) 


The  expected  values  are  obtained  by  multiplying  the  quantities  in 
square  brackets  by  the  joint  prior  density  P(2i»  •  •  •2inl™'^*'”^^ 
given  in  (42)  and  integrating  over  the  validation  gate. 


A 

The  second  expression  in  (21)  is  obtained  as  a  intermediate 
result  in  Appendix  B. 
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Considerable  simplifications  result  if  one  applies  a  linear 
transformation  of  variables  to  obtain  a  spherical  gate, 
followed  by  a  change  to  spherical  coordinates.  Because  of  the 
spherical  symmetry  of  the  gaussian  density  and  of  the  expressions 
(20)-(21),  off-diagonal  elements,  cross-terms,  and  angular 
variables  drop  out  like  flies,  leaving  scalar  integrals  over  the 
radial  variables  (i.e.,  over  Il2jll»  j=l,...m).  The  detailed 
derivation,  which  is  given  in  Appendix  B,  leads  to 


t 

4 

< 

f 

i 

i 


i 


Uj^(m)  =  m 


U2(m)  =  m 


M 


PpP^m  +  (l-PnPrt)CV  (27r)M/2 


iii 


PpPgm  +  (1-Pj3Pq)CV  (2 


I2  (n»)S 


(22) 


(23) 


where  the  scalar  integrals  and  l2(in)  are  defined  as 

g 

-  J"  r**'’’^exp(-r^/2)dr 


0 

g 


g  exp{-rJ)rj^2 

12C)  «  — - - 

0  0  b  +  ^exp(-r  j2/2) 


(24) 


(rj_r2. .  ^dr]_. 


.dr 


m 

(25) 


and  the  constant  b  is  defined  in  (13) . 

Substituting  (22) -(25)  and  (37)  into  (18)  and  cancelling 
leads  to 


£k  •  (qi-q?)wsw" 


(26) 
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where 


(2  7t)M/2  (n,_l)l  (2tt)”/2 


24 _  ^ 

.  .m/2  ^ 


00  e"^^(CV)"*”^ 


It  =  Pr 


Cm 


(27) 


Cj4  «  e"CV(CV)"*"^  /  mV""^ 

‘32  =  Pd  ---M/i  E  -  1-m)  l2("') 


(28) 


and  substitution  of  this  and  (17)  into  (14)  yields 


£k|k  -  £k|k-i  -  (wsw')  (29) 

This  can  be  simplified  further  by  noting  that  for  typical 
values  of  g  and  M  (g=4  or  5  and  M<10) ,  Pq  is  approximately  1  and 
is  approximately  Pq  (substitute  x^-/^  for  r  in  to  get  a 
gamma  function) . 

The  upshot  of  all  this  is  that  the  deterministic 
approximation  to  the  covariance  equation  becomes 


£k|k  -  £k|k-l  -  92<«SW') 

where  the  scalar  parameter  q2»  defined  in  (28)  ,  depends  on  Pq  and 
C  and  lies  between  0  and  1.  Comparing  this  to  (8) ,  it  is  clear 
that  the  factor  q2  reduces  the  covariance  improvement  due  to  the 
term  WSW':  the  smaller  q2  is,  the  greater  the  degradation. 
Thus,  q2  provides  a  useful  measure  of  tracking  performance  as  a 
function  of  Pq  and  C. 


Evaluation  of  q2 

A  Monte  Carlo  integration  program  has  been  used  to  integrate 
(25)  and  compute  q2  for  various  values  of  Pjj  and  C.  Evaluating 
the  integrand  in  (25)  at  5000  random  points  and  truncating  the 
sum  in  (28)  at  10  terms  yields  a  table  of  q2  values  with  error 
standard  deviations  of  0.5-1. 5%  for  M*1  and  2,  and  1-5%  for  M=3. 

It  is  useful  to  display  these  numbers  via  contour  plots  of 
q2f  as  shown  in  Figures  2-4  (the  curves  become  smoother  if  more 
points  are  used  in  the  Monte  Carlo  integration)  .  The  vertical 
axis  is  Pq  and  the  horizontal  axis  is  CV,  the  expected  number  of 
false  measurements  in  the  gate  (proportional  to  false  alarm 
probability) .  A  4-sigma  gate  (g»4)  was  used  and  the  three  plots 
correspond  to  measurement  vectors  ^  with  dimensions  M  »  1,  2,  and 
3,  respectively.  We  shall  refer  to  these  as  tracker  operating 
characteristic  (TOC)  curves,  by  analogy  with  the  receiver 
operating  characteristic  (ROC)  curves  to  be  given  in  Section  4 
below. 


Note  that  V  =  §.1  £  ~  ^^]c|k  that  the 

gate  volume  in  fact  depends  upon  the  state  error  covariance  matrix 
— kIk-1*  follows  that  the  above  analysis  is  fully  valid  only 

in  steady  state  ~  constant).  We  are  currently  working  on 

extending  this  procedure  by  iterating  equations  (4),  (30),  and 
(28)  to  obtain  P(Pjj,C)  ,  the  steady-state  error  covariance  matrix. 

A  slightly  modified  set  of  TOC  curves  can  then  be  formed  with 
|p1  '  (or  some  other  appropriate  metric)  plotted  as  a  function 
of  Pp  and  C.  In  the  process  of  obtaining  the  table  of  P(Pjj,C)  , 
we  also  expect  to  determine  empirically  those  values  of  Pj^  and 
C  for  which  the  covariance  equations  become  unstable. 
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4.  ANALYSIS  OP  DETECTION  ALGORITHMS 

The  tracking  inputs  shown  in  Figure  1  can  come  from  a 
variety  of  sources,  ranging  from  simple  threshold  detectors  to 
complex,  multi-dimensional  signal  processing  algorithms. 
Nevertheless,  it  is  often  possible  to  view  the  source  of 
measurement  data  as  a  classic  receiver  structure  in  which  a 
threshold  level  and/or  other  parameters  control  the  probabilities 
of  detection  (Pq)  and  false  alarm  (Pp^) • 

In  passive  sonar  applications,  for  example,  time  series  data 
from  each  sensor  is  often  passed  through  a  beamformer,  resulting 
in  several  data  streams  (beams)  corresponding  to  different  angles 
of  arrival  (azimuths)  ,  and  then  a  sequence  of  fast  Fourier 
transforms  (FFTs)  is  performed  on  each  data  stream.  The  FFT 
operation  may  be  viewed  as  a  bank  of  filters,  so  that  one  deals 
with  many  separate  narrowband  time  series,  one  from  each  beam-bin 
or  resolution  "cell"  in  the  two-dimensional  frequency-azimuth 
space . 


There  are  many  different  detection  algorithms  for  dealing 
with  these  data  streams;  their  structures  depend  upon  the  nature 
of  the  application  and  the  expected  signals.  A  simple  but  common 
procedure,  if  the  signals  are  narrowband  and  relatively  stable, 
is  to  form  power  spectra  (i.e.,  take  the  squared  magnitude  of 
each  sample)  and  then  integrate  (i.e.,  sum)  each  data  stream  over 
an  appropriate  period  of  time.  At  the  end  of  each  integration 
period,  a  noise  spectrum  equalization  (NSE)  procedure  compensates 
for  different  background  noise  levels  at  different  frequencies, 
and  detections  are  declared  in  all  beam-bins  for  which  the 
integrated  output  exceeds  a  threshold. 

Except  for  the  NSE,  this  amounts  to  passing  each  data  stream 
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through  an  energy  de tec tor  [19,  20]  ,  which  assumes  an  unknown 
signal  in  white,  gaussian  noise.  Under  the  noise-only 

hypothesis,  the  test  statistic  has  a  chi-square  distribution  with 
2TW  degrees  of  freedom,  where  T  is  the  integration  time  and  W  is 
the  bin  width.  With  the  signal  present,  it  has  a  non-central 
chi-square  distribution  with  the  same  degrees  of  freedom  and 
non-centrality  parameter  equal  to  the  SNR. 

We  shall  use  these  distributions  to  compute  receiver 
operating  characteristics  (ROCs) .  A  standard  ROC  curve  is  the 
locus  of  points  obtained  as  the  detection  threshold 

varies  over  all  possible  settings,  and  a  family  of  such  curves 
for  different  SNRs  is  usually  plotted. 

In  this  case  an  energy  detector  is  operating  on  each 
frequency-azimuth  cell,  and  we  are  interested  in  the  expected 
number  of  false  alarms  in  all  cells  of  interest,  which  is  the 
single-cell  false  alarm  probability  multiplied  by  the  number  of 
independent  cells  in  a  tracker's  validation  gate.  Most 
algorithms  require  detections  to  be  separated  by  one  or  more 
empty  cells  (i.e.,  adjacent  cells  above  threshold  are  considered 
part  of  the  same  detection)  ,  and  the  effective  number  of  cells 
generally  depends  in  a  complicated  way  on  the  details  of  the 
algorithm.  For  the  sake  of  this  discussion,  we  shall  simply 
assume  that  false  alarms  can  potentially  occur  in  1/3  of  the 
cells  in  any  given  integration  period. 

A  typical  automatic  line  detection  algorithm  of  the  type 
described  above  might  have  an  integration  time  of  250  sec.  and 
FPT  resolution  of  0.1  Hz,  or  50  degrees  of  freedom,  and  each  cell 
would  be  1  beam  by  0.1  Hz.  The  tracker  in  this  situation  might 
have  an  ' innovation  covariance  matrix 


S 


1  beam 
0 


0.3  HzJ 


(31) 
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and  a  4-sigma  validation  gate  with  volume 


V  =  C2g^|S|^/2  ^  t;(16)  (0.3)  =  15  beam-Hz  (32) 

The  total  number  of  resolution  cells  in  the  gate  is  thus 
15/0.1  =  150,  and  the  number  of  independent  cells  is  150/3  =  50. 
It  follows  that  the  expected  number  of  false  alarms  per  gate  is 
CV  =  SOPjpj^. 

A  set  of  ROC  curves  for  an  energy  detector,  computed  using 
the  above  assumptions,  is  shown  in  Figure  5.  For  a  given  SNR, 
one  can  adjust  the  detection  threshold  so  as  to  select  an 
operating  point  anywhere  along  the  corresponding  ROC  curve. 


FIG.  5.  TYPICAL  ROC  CURVES  FOR  ENERGY  DETECTOR 

More  sophisticated  narrowband  line  detectors  and  line 
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trackers  utilize  complex  spectra,  and  their  effective  integration 
times  are  usually  less  than  those  described  above.  Although  the 
details  of  these  algorithms  vary  greatly,  their  basic  operations 
tend  to  be  similar  to  those  of  a  quadrature  receiver  or 
incoherent  matched  filter  [20] ,  which  assumes  a  sinusoidal  signal 
of  unknown  phase.  Under  the  signal-plus-noise  hypothesis,  the 
test  statistic  has  a  Rician  distribution,  which  reduces  to  a 
Rayleigh  distribution  in  the  noise-only  case. 

A  set  of  ROC  curves  for  a  quadrature  receiver  is  shown  in 
Figure  6.  These  were  computed  using  the  expressions  given  in 

[20]  for  Pp^  and  Pj^,  and  using  the  same  validation  gate  as 
above.  Another  set,  for  the  case  where  both  amplitude  and  phase 
are  unknown  (slow  Rayleigh  fading),  is  shown  in  Figure  7. 


FIG.  6.  ROC  CURVES  FOR  QUADRATURE  RECEIVER 

Another  class  of  algorithms  for  passive  sonar  applications 
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18  dB 


FIG.  7.  ROC  CURVES  FOR  QUADRATURE  RECEIVER  WITH  RAYLEIGH  FADING 

( 

involves  coherent  or  incoherent  correlation  of  data  received  at 
two  geographically  distinct  sites.  The  resulting  statistic,  as  a 
function  of  time  delay  and  frequency  (Doppler)  shift,  is  known  as 
an  ambiguity  surface ,  and  a  signal  received  at  both  sites  from 
the  same  target  produces  a  peak  in  the  surface.  Detection  of 
peaks  above  some  threshold  is  qualitatively  similar  to  the 
situation  described  above,  in  which  peaks  are  sought  in  the 
frequency-azimuth  plane.  With  appropriate  assumptions, 

approximate  ROC  curves  can  be  obtained  for  such  receivers. 
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5.  OPTIMIZATION  OF  TRACKING  PERFORMANCE 

Superimposing  one  of  the  sets  of  TOC  contours  from  Figures 
2-4  on  one  of  the  ROC  curves  from  Figures  5-7,  an  operating  point 
(i.e.,  threshold)  can  be  selected  so  as  to  maximize  q2  and 
"optimize"  tracking  performance.  For  example.  Figure  8 
shows  the  TOC  contours  for  M=2  superimposed  on  the  ROC  curve 
corresponding  to  a  quadrature  receiver  and  SNR  =  8  dB.  The 
"optimal"  point  is  indicated  by  0. 


Cluittr  tfmwiy  CV  -  50 


FIG.  8.  OPTIMAL  THRESHOLD  SETTING  FOR  QUADRATURE  RECEIVER  WITH 
SNR  =  8  DB 

This  graphical  method  for  selecting  an  operating  point  can, 
of  course,  be  replaced  by  a  mathematical  optimization;  the 
obvious  necessary  condition  is  that  the  ROC  and  TOC  curves  be 
tangent.  However,  the  practical  difficulty  of  computing  the 
required  differentials  to  solve  the  necessary  conditions  is 
substantial . 
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Another  interesting  issue  is  the  optimization  of  tracking 
performance  when  the  signal's  SNR  is  not  known.  In  this  case, 
several  ROC  curves  are  involved  and  one  must  select  a  threshold 
that  is  best  (in  some  sense)  for  a  whole  range  of  SNRs. 
Alternatively,  an  adaptive  scheme  could  be  devised,  whereby  the 
SNR  is  monitored  and  the  threshold  adjusted  so  as  to  maximize  q2 
along  the  current  ROC  curve. 
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6.  APPLICATION  TO  MODIFIED  LIKELIHOOD  FUNCTION 

Calculation  of  the  modified  likelihood  function  [18] 

requires  the  second  moment  U2(m)  of  the  combined  innovation  (10) , 

k—  1 

conditioned  on  Y  and  m.  This  may  be  written  as 

U2(m)  =  q2{ni)S  (33) 

where  q2(n>)  represents  all  the  scalar  factors  on  the  right-hand 
side  of  (23) . 

Since  this  is  proportional  to  the  single-measurement 
innovation  covariance  matrix  S,  the  standard  likelihood  function 
machinery  need  only  incorporate  the  scale  factor  q2(ni)  in  order 
to  compute  the  modified  version.  Values  of  q2(®)  various  m, 

g,  Pjj,  and  C  may  be  precomputed,  using  the  methods  described  in 
Section  3  above. 
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7 .  CONCLUSION 


We  have  explored  an  important  connection  between  data 
association/tracking  algorithms  and  the  upstream  signal 
processing/detection  algorithms  that  provide  their  measurements. 
A  quantitative  relationship  has  been  derived  that  specifies  the 
effect  on  tracking  performance  of  threshold-dependent  parameters 
such  as  the  probabilities  of  detection  and  false  alarm,  and  a 
graphical  technique  has  been  described  for  selecting  an  optimal 
threshold. 
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APPENDIX  A 

PROBABILITY  CALCULATIONS 


In  this  appendix  we  shall  derive  a  number  of  expressions 
needed  in  the  main  text.  Letting  Yj(m)  denote  the  prior 
probability  of  the  event  Xj»  conditioned  on  m,  the  total 
probability  theorem  yields 

Yj  (m)  »  P{xj  |m,Y'«“l}  »  P{xj|m} 

=  P{x  j  ljn^»m-l,m}p{m^»m-l|m}  +p{xj  l®^*tn,m}p{m^»mlm} 

I  (1/m) P{m^«m-1 |m}  +  (0) P{m^«m|m} ,  j»l,..,m 

“  I  f  P  >1  r  P  .  1 

'  (0)  P{m^*m-1  |m}  +  (1)  p{m^*m |m} ,  j»»0 

because  m^ ,  the  number  of  false  measurements,  must  be  either  m-1 
(if  the  target  is  detected)  or  m  (if  it  is  not)  .  Using  Bayes' 
rule  and  the  assumed  Poisson  distribution  for  false  measurements, 

p{m^=m-l|m}  =  p{m|m^»m-l}p{m^»m-l}/p{m} 

»  [PijPq]  [e-^(CV)"'"V(m-l)  I]/P{m} 

-  Pl3PGm/[PuPQm  +  (l-PpPQ)CV]  (35) 

P{m^*m|m}  =  p{m|m^»m}p{m^am}/p{m} 

-  [1-PdPq]  (e-CV(cv)m/n,ij/p|n,| 
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(l-PDPG)CV/tPuPGm  +  (1-Pj3Pq)CV] 


(36) 


where  the  denominator  P{m}  is  the  prior  probability  of  m  and  is 
equal  to  the  sum  of  the  numerators  in  the  two  equations: 


P{m}  =  p{m(Y'«“^} 

=  [P^Pcm  +  (l-PDPG)CV]e-CV(cv)"'"Vm!  ,  m*0,l,...  (37) 


Substituting  back  into  (34)  yields 


Yj  (m)  « 


I  +  (I-PdPg)CV], 

I (l-PDPQ)CV/rPoPQm  +  (I-PuPq)CV], 


3®lf • • •m 

j»0 


(38) 


Note  that  Yj (m)  is  independent  of  j  for  j>0. 

Using  Bayes'  rule,  the  posterior  probabilities  in  (10)  can 
be  expressed  as 


6j  £  P{Xj|Y''}  =  P{Xjl2i....2m'"‘'Y'''^} 


=  P(2i, 

“  P(2ir 


•2inl  Xj  fin,Y''"l 


•  inilx  j  rn'rY''"! 


)P{Xj|in,Y'^“l}/p(2l,.. 

numerators 


.2^|m,Y'^-M 


(39) 


The  first  factor  is  the  joint  probability  density  of  the  m 
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candidate  measurements,  conditioned  on  the  j-th  one  being 
correct.  According  to  the  PDA  assumptions,  the  correct 
measurement  £j  has  a  gaussian  density 

N(2j;0,S)/PQ  =  (l/PQ)exp(-£1S‘^2j/2)/(2TT)”/2|s|l/2  ^40) 

with  mean  £  and  covariance  S,  where  the  factor  1/Pq  accounts  for 
its  restriction  to  the  validation  gate,  and  each  incorrect 
measurement  has  a  uniform  density  V“^.  It  follows  that 


I  V-m+lN(2j;0,S)/PQ,  j  =  l 

j=0 


.m 

(41) 


The  second  factor  in  (39)  is  the  prior  probability  of  Xj,  given 
by  (38)  .  The  denominator  is  the  joint  prior  density  of  the 
measurements,  conditioned  only  on  m  (and  the  past  data), 


=  (42) 

m 

V^'Yodn)  +  V~™+^^(l/PQ)N(2j;0,S)  Yj  (m) 

Note  that  with  the  above  conditioning,  the  validated  measurements 
are  not  independent,  i.e.,  (42)  is  not  equal  to  the  product  over 
j  of  the  marginal  prior  densities 

p(2j  |m,Y'^-l)  =  V”^[l-Yj(m)]  +  (1/Pq)  N  (2j  ?0  ,S)  Yj  (m)  (43) 

Finally,  substitution  of  (38)  ,  (41)  ,  and  (42)  into  (39) 

followed  by  a  certain  amount  of  rearrangement  yields 
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j=l, . . . 


(44) 


where 

b  S  (2Tr)M/2c|s|l/2(i_p^p^)/p^ 

=  (27r)M/2(cv/c„gM)  (1-PdPg)/Pd 


in 


(45) 


(46) 
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appendix  b 

DERIVATION  OF  Dj^  AND  O2  INTEGRALS 


In  this  appendix  the  expectations 
Uj^(m)  ^  e.y.y' 


(B-1) 


and 

U2(m)  ^  ^i^i  .1^  (B-2) 

will  be  evaluated,  where  is  given  by  (44)  and  the  joint  density 
of  is  (42).  In  order  to  simplify  the  arguments,  we 

shall  make  use  of  the  fact  that  S  is  positive  definite  and  change 
variables  to 


i  ~  l,«..m 


(B-3) 


so  that 


y's“^y 


i  ='^i'^i 


(B-4) 


and 


^>2 


y  .y'  .  s'^S.t.S 


(B-5) 
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In  terms  of  the  new  variables  the  validation  gate  becomes 

a  sphere  :  lU  H  ^  g^)  with  volume  c^^g^  =  V;  the  B^^'s  can 

be  rewritten  as 

A  ^ 

8.  =  - - - 5—  ,  i  =  l,...m  (B-6) 

^  r  -%lk  Ir 

b  +  I  e  ^  j  '' 

j=l 


where 


b^  (27r)”/^(CV/Cj^g”)  (1-PjjPq)/Pj^ 


(B-7) 


is  the  same  as  in  (46) .  Note  that  the  dimensionless  quantity  CV, 
the  average  number  of  false  measurements  in  the  (now  spherical) 
gate,  is  unaffected  by  the  variable  change.  The  joint  density 
(42)  becomes 


P(^ 


1' 


m 


m,Y  )  = 


Yq 


-m  ^  ,  .,''-m+l 

+  Yj(ni)  V 


m 

I 

i=l 


;/(!;  .  ;0,I)/Pp^  (B-8) 


Ikill  1  ^  ®  otherwise. 
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Using  the  change  of  coordinates  from  y  to  C ,  expressions  (19) 
through  (21)  may  be  reexpressed  as  follows: 


yk-l] 

=  WS^[U^(m)-G2(ni)]s’^W'  ,  m=l,2,... 

(B-9) 

(m) 

=  E[  ^  !;(|z^"-^,m]  =  S“'^U  (m)s"^ 

i=l  1  1  1  ^ 

(B-10) 

U2(m) 

in  ^  A  T  j  1 

=  E[  1  &  ^  1  6  !z^"^,m]  =  s"""U2(m)s"^ 

i=l  ^  ^  j=l  J  ^ 

(B-11) 

Let  D  stand  for  the  spherical  validation  gate  region.  Then  the 
above  two  expectations  may  be  written  out  explicitly  as  multiple 
integrals  of  certain  matrices,  taken  over  m  copies  of  D; 


^1  = 


^2  = 


* 

•  •  •  « 

J 

J 

D  D 


A(^2^,  .  *  *^^m  (B"”! 2 ) 


B(C, .,^„)P ,c„)d;  . . .d;^  (B-13) 

1  mi  m  1  m 


D  D 


where  p  is  the  joint  innovations  density  (B-8)  with  the  conditioning 
suppressed  and  A,B  are  Mxm  matrices  described  componentwise  by 


m 


pq 


iil  ^ 


1  £  p,q  £  M 


(B-14) 


m 


B 


pq  - 


1  <  p,q  <  M 


(B-15) 
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Here  the  M-dimensional  vector  has  been  denoted  as 
The  integrals  (B-12)  and  (B-13)  are  actually  mM-fold  integrals. 
Although  the  complexity  of  undertaking  such  a  large  computation 
directly  would  prove  formidable,  a  number  of  observations  show 
that  these  integrals  have  a  fairly  simple  and  straightforward 
structure: 


Observation  1.  The  matrices  and  U2  are  diagonal. 

From  the  expressions  (B-14)  and  (B-15)  one  may  deduce  that 

off-diagonal  elements  will  integrate  out  to  zero.  For  when 

p  ^  q,  both  Apq  and  Bp^  become  odd  polynomials  of  second  degree 

in  the  variables  coefficients  either  single 

6's  or  pairwise  products  of  S's.  Now  both  the  B's  and  the 

two  terms  of  the  probability  density  p  in  (B-8)  are  positive 

functions  depending  only  on  /  •  •  •  41  ,  i.e.  they  have 

the  scune  values  at  antipodal  points  of  a  sphere.  The  odd 

polynomials  from  A  _  and  B  _  will  have  opposite  values  at 

pq  pq 

antipodal  points,  so  that  their  contributions  to  the  total 
integral  will  cancel. 


Physically/  this  amounts  to  the  observation  that  off  -diagonal 
elements  of  the  inertia  tensor  vanish  in  a  principal  axis  system, 
for  a  spherical  body  with  the  shape  of  the  validation  region  and 
with  mass  distribution  density  appropriately  defined  in  terms  of 
the  0's  and  p. 
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Observation  2.  Integrals  of  cross  terms  (i^j)  vanish  in  the 
svim  (B-15)  defining  B 

pq 

From  the  previous  observation,  we  need  only  concern 
ourselves  with  diagonal  entries  Bpp*  Again  when  ij^j  we  have  an 
odd  polynomial  in  the  ^'s,  and  the  same  arguments  as  before 
apply  again  to  show  its  integral  over  a  sphere,  weighted  by 
the  density  p,  must  vanish. 

From  the  point  of  view  of  the  original  definition  (B-11) 
of  U2f  we  could  say  that  distinct  innovations  are  orthogonal; 

their  inner  products  vanish  and  do  not  contribute  to  the  second 
moment  of  the  combined  innovations. 


Observation  3.  Each  term  of  A  or  B  making  a  nonvanishing 
contribution  to  N=1  or  2,  has  the  same  value,  given  explicitly 
by 


a 

r 


o  r 


=  K 


j=l 


dr 


m 


(B-16) 


where 


=  C„(m,M)  ^  ^  -■  y 

^  ^  (2TT)”/2(c„g”)"'"^P, 


Yj  (m) 

57r 


(B-17) 
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(B-18) 


C  =  C  (m,M)  ^  ^ 


K  =  K(m,M)  = 


(B-19) 


To  prove  this,  let  us  consider  the  first  term  of  the  (1,1) 

12 

entry  of  each  integral,  which  cam  be  written  as  .  ,?jn)  , 

N=l,2.  Since  p  and  depend  only  on  the  norms  [|  ^  j^||  ~ 

it  is  clearly  advantageous  to  introduce  spherical  coordinates 
in  each  of  the  m  copies  of  the  ^-space: 


d?.  *  d^Jd;^...d!;”  *  r^’^dr^dco^"^ 


(B-20) 


Here  dcu^”^  denotes  the  "unit  solid  angle,"  i.e.  the  surface 


1  ^  o  o  M  ^ 

element  on  the  unit  sphere  +  ...  +  (c")  =  1. 


Making  these  substitutions,  the  integral  for  the  expected 


"'N  12 

value  of  becomes 


e 


-hri 


2  nN 


?  -Jsr? 
b+  e 


j=l 


M-1, 


M-1, 


,M-1, 

'm 


M- 


{B-21) 
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The  constants  Cg,C^,  and  K  defined  in  (B-17) ,  (B-18)  are 
introduced  to  simplify  the  notation. 


The  expression  in  (B-21)  is  not  yet  complete  because  the 
variable  still  has  to  be  reexpressed  in  angular  co^^-coordinates 
for  the  first  copy  of  the  M-dimensional  state  space.  However, 
the  other  angular  variables  occur  in  the 

integrand,  so  that  their  integrals  can  be  absorbed  in  the 
constant  K.  The  fact  needed  to  do  this  is  that  the  integral 
of  the  unit  solid  angle  is  just  the  surface  area  of  the  unit 
sphere,  which  is  M  times  c„  (the  volume  of  the  unit  sphere) . 


unit  sphere 


To  handle  the  C ^-dependence 
of  the  integrand  in  (B-21)  we  are 
going  to  use  an  argument  almost 
identical  to  the  one  used  to  derive 
the  volume  of  the  unit  sphere 
(which  is  to  start  from  a  known 
integral  and  "work  backwards"). 

The  known  integral  that  we  wish  to 
exploit  is 


(cj;)^exp(-Js  I  C^)dcJdC?...dc”  = 


\Z< 


(B-22) 


Now  if  this  integral  is  rewritten  in  spherical  coordinates,  we 
should  make  a  substitution  =  r.^f  where  f  is  a  smooth 

function  independent  of  r^^.  For  example  if  coordinates  are 
chosen  on  the  sphere  as  depicted  in  the  figure,  we  would  have 
f(a)^”^)  =  cos  0,  where  the  angle  6  is  measured  downward  from  the 
north  pole  as  shown.  Rewriting  (B-22)  in  spherical  coordinates 


gives 


(2,r)«/2  ^  ^M-1  ,2  ^^M-1 


(B-23) 


sphere 


The  substitution  Hr^  =  transforms  the  radial  part  of 
(B-22)  into  the  well-known  integral  representation  for  the 
gamma  function 


r(t)  =  e"*l  xj"^dx^ 


(B-24) 


so  now  we  have 


I,  I  £2, <,«-!) 


(B-25) 


sphere 


^7e  are  now  in  a  position  to  solve  for  the  nonradial  part  of  the 
integral  in  (B-25) .  But  when  this  is  done,  the  quantity  we 
obtain  is  recognized  to  be  c^^,  the  volume  of  the  unit  sphere: 


.2,  M-1,,  M-1 


TT 


M/2 


r(|  + 


=  c, 


1) 


M 


(B-26) 


Actually,  by  using  properties  of  the  gaunma  function,  this 
expression  for  c^^  can  be  written  more  explicitly  as 


-M/2  M/2  M/2 

2  TT  _  IT 

— -  -  -  if  M  IS  even 

M(M-2)...4*2  (m/2)! 


'M 


il 

M  ^M-2 


2(M+1)/2  ^(M-1)/2 
M(M-2) . . .5»3«1 


if  M  is  odd 


(B-27) 


Now  returning  to  the  original  integral  (B-21) ,  we  note  that 
even  though  the  integrand  is  not  the  same  as  {B-22) ,  the  angular 
dependence  so  that  (B-26)  may  be  applied.  Absorbing  this 
c,,  into  K  shows  that  the  6^  (c][)  ^-contribution  to  is  given 
exactly  by  (B-16) .  It  should  be  noted  that  this  formula  is 
also  immediately  valid  for  M*1  without  even  introducing  the 
spherical  change  of  coordinates  (B-20) ,  since  the  volume  of  the 
unit  interval  is  2  (correctly  given  by  B-27) . 


Although  we  began  the  discussion  by  considering  the  contri- 

1  2 

bution  of  a  single  term  to  the  expectation,  the 
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homogeneity  of  the  expression  (B-21)  shows  that  the  other  terms 

^  in  the  (1,1) -entry  will  each  make  an 
identical  contribution.  The  same  argument  applies  to  other 
diagonal  entries  of  and  U2S  homogeneity  again  shows  that 
each  term  will  make  the  same  contribution,  given  by  (B-16) .  This 
justifies  the  conclusion  that 

^  (B-28) 

where  I  is  the  identity  matrix. 

Next,  let  us  transform  back  from  ^ -coordinates  to  the 
original  v-coordinates  as  given  in  (B-5) .  We  have  proven  that 
our  original  expectations  Uj^,U2  are  scalar  multiples  of  the 
covariance  matrix  S: 

Ujj  =  S^(mnjjI)S*®  =  mrijjS  (B-29) 

Finally,  substitution  of  (B-17) - (B-19)  and  (38)  into 
(B-16)  followed  by  numerous  cancellations  yields  the  expressions 
(22)-(25)  in  Section  3. 
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